

*************************************************************************************************
* Most common specialty combinations 
*************************************************************************************************
* Physician - Emergency 
use "${Gpath}\data\\proceduralist_physician_level.dta",  clear 

keep code1 code2 cod_olap  sp_naih tot_vi
ren sp_naih n_aih 
merge m:1 n_aih  using "${Gpath}\data\\estimation_sample.dta", nogen keepusing(n_aih main_sample) keep(3)
keep if main_sample==1 
gen cases=1 
set seed 123 
gen random=runiform()
gsort n_aih -tot_vi random
by n_aih: keep if _n==1 
collapse (sum) cases  , by(code1) 
sum cases
gen share=cases/r(sum)
keep code1   share 
gsort -share 
keep if _n<=10 
gen id=_n 
ren share share_physician 
replace share_physician=100*share_physician
replace code1="Cardiology + Internal Medicine" if _n==1 
replace code1="Cardiology" if _n==2 
replace code1="Cardiology + Internal Medicine + Critical Care Medicine" if _n==3
replace code1="Internal Medicine" if _n==4
replace code1="Cardiology + Critical Care Medicine" if _n==5
replace code1="Cardiology + Internal Medicine + General Surgery" if _n==6
replace code1="Cardiology + Internal Medicine + Family Medicine" if _n==7 
replace code1="Cardiology + Internal Medicine + Cardiovascular Surgery" if _n==8
replace code1="Internal Medicine + Preventive Medicine and Public health" if _n==9 
replace code1="Cardiology + Internal Medicine+ Family Medicine + Preventive Medicine and Public health" if _n==10 
ren code1 specialty_physician 
save A.dta, replace 

* Proceduralist - Emergency 
use "${Gpath}\data\\proceduralist_physician_level.dta",  clear 
keep code1 code2 cod_olap  sp_naih tot_vi
ren sp_naih n_aih 
bysort n_aih: keep if _n==1 
merge m:1 n_aih  using "${Gpath}\data\\estimation_sample.dta", nogen keepusing(n_aih main_sample) keep(3)
keep if main_sample==1 
gen cases=1 
gsort n_aih -tot_vi
by n_aih: keep if _n==1 
collapse (sum) cases , by(code2) 
sum cases
gen share=cases/r(sum)
keep code2   share 
gsort -share 
keep if _n<=10 
gen id=_n 
ren share share_proceduralist 
replace share_proceduralist=100*share_proceduralist

replace code2="Cardiology + Internal Medicine" if _n==1 
replace code2="Cardiology + Internal Medicine+ Radiology" if _n==2 
replace code2="Cardiology" if _n==3 
replace code2="Cardiology + Radiology" if _n==4
replace code2="Cardiology + Internal Medicine+ Cardiovascular Surgery" if _n==5
replace code2="Cardiology + Internal Medicine+ Anesthesiology" if _n==6
replace code2="Cardiology + Cardiovascular Surgery" if _n==7 
replace code2="Cardiology + Internal Medicine+ Critical Care Medicine" if _n==8 
replace code2="Cardiology + Internal Medicine+ Anesthesiology + Critical Care Medicine" if _n==9
replace code2="Cardiology + Internal Medicine+ Anesthesiology + Critical Care Medicine + Radiology" if _n==10  

ren code2 specialty_proc 
merge 1:1 id using A.dta, nogen keep(1 3)
erase A.dta 
order specialty_physician share_physician 
drop id 

outsheet using "${Gpath}\table\\common_specialties.csv", replace  

* Physician - Non-Emergency 
use "${Gpath}\data\\proceduralist_physician_level.dta",  clear 
keep code1 code2 cod_olap  sp_naih tot_vi
ren sp_naih n_aih 
merge m:1 n_aih  using "${Gpath}\data\\estimation_sample.dta", nogen keepusing(n_aih main_sample) keep(3)
keep if main_sample==0
gen cases=1 
set seed 123 
gen random=runiform()
gsort n_aih -tot_vi random
by n_aih: keep if _n==1 
collapse (sum) cases  , by(code1) 
sum cases
gen share=cases/r(sum)
keep code1   share 
gsort -share 
keep if _n<=10 
gen id=_n 
ren share share_physician 
replace share_physician=100*share_physician

replace code1="Cardiology + Internal Medicine" if _n==1 
replace code1="Cardiology" if _n==2 
replace code1="Cardiology + Internal Medicine + Critical Care Medicine" if _n==3
replace code1="Cardiology + Internal Medicine + Family Medicine" if _n==4 
replace code1="Cardiology + Internal Medicine + General Surgery" if _n==5 
replace code1="Cardiology + Critical Care Medicine" if _n==6 
replace code1="Internal Medicine" if _n==7 
replace code1="Cardiology + Internal Medicine + Cardiovascular Surgery" if _n==8 
replace code1="Cardiology + General Surgery" if _n==9 
replace code1="Internal Medicine+General Surgery+ Family Medicine + Preventive Medicine and Public health" if _n==10  
ren code1 specialty_physician 
save A.dta, replace 

* Proceduralist - Non-Emergency 
use "${Gpath}\data\\proceduralist_physician_level.dta",  clear 
keep code1 code2 cod_olap  sp_naih tot_vi
ren sp_naih n_aih 
bysort n_aih: keep if _n==1 
merge m:1 n_aih  using "${Gpath}\data\\estimation_sample.dta", nogen keepusing(n_aih main_sample) keep(3)
keep if main_sample==0 
gen cases=1 
gsort n_aih -tot_vi
by n_aih: keep if _n==1 
collapse (sum) cases , by(code2) 
sum cases
gen share=cases/r(sum)
keep code2   share 
gsort -share 
keep if _n<=10 
gen id=_n 
ren share share_proceduralist 
replace share_proceduralist=100*share_proceduralist

replace code2="Cardiology + Internal Medicine" if _n==1
replace code2="Cardiology + Radiology" if _n==2
replace code2="Cardiology + Internal Medicine+ Radiology" if _n==3
replace code2="Cardiology" if _n==4
replace code2="Cardiology + Cardiovascular Surgery" if _n==5
replace code2="Cardiology + Internal Medicine+ Cardiovascular Surgery" if _n==6
replace code2="Cardiology + Internal Medicine+ Anesthesiology" if _n==7
replace code2="Cardiology + Critical Care Medicine + Cardiovascular Surgery" if _n==8
replace code2="Cardiology + Critical Care Medicine" if _n==9
replace code2="Cardiology + Internal Medicine+ Critical Care Medicine" if _n==10

ren code2 specialty_proc 
merge 1:1 id using A.dta, nogen keep(1 3)
erase A.dta 
order specialty_physician share_physician 
drop id 

outsheet using "${Gpath}\table\\common_specialties_nonemergency.csv", replace  

************************************************************************************************* */
* Baseline estimation
************************************************************************************************* */
use  "${Gpath}\data\\estimation_sample.dta", clear 
keep if main_sample==1 

estimates drop _all 
local h=0 

qui sum m_30d
scalar Mean=r(mean)

* Baseline specification 
qui reghdfe  m_30d  share_dcbo_* overlap, a(cpfxcnesxyear) cluster(sp_cnes)
eststo M`h++', addscalars(Mean Mean )

* hospital x (month- and day-of-week) FE 
qui reghdfe  m_30d  share_dcbo_* overlap , a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsing
eststo M`h++', addscalars(Mean Mean )

* patient characteristics
qui reghdfe  m_30d  share_dcbo_* $XVar0 overlap , a(cpfxcnesxyear) cluster(sp_cnes) 
eststo M`h++', addscalars(Mean Mean )

* Physician characteristics
qui reghdfe  m_30d  share_dcbo_* $Xvar0_phy overlap, a(cpfxcnesxyear) cluster(sp_cnes) 
eststo M`h++', addscalars(Mean Mean )

* All 
qui reghdfe  m_30d  share_dcbo_* $XVar0 $Xvar0_phy overlap, a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsing
eststo M`h++', addscalars(Mean Mean )

esttab M* using "${Gpath}\\table\\mortality_baseline.csv",   keep( *overlap) replace star(* 0.1 ** 0.05 *** 0.01) ///
se b(4) staraux brackets   ///
compress   begin(";")  delimiter(";") end(";")  scalar(Mean ar2 r2_a N )



************************************************************************************************* */
* Controlling for specialty and adjusting the specific overlap measure 
************************************************************************************************* */
use  "${Gpath}\data\\estimation_sample.dta", clear 
keep if main_sample==1 

estimates drop _all 
local h=0 
qui sum m_30d
scalar Mean=r(mean)

qui reghdfe  m_30d  overlap  share_dcbo_* $XVar0 $Xvar0_phy , a(cpfxcnesxyear  hospitalxdow hospitalxmonth)   resid cluster(cnes)
eststo M`h++', addscalars(Mean Mean )

qui reghdfe  m_30d  overlap  rw*  $XVar0 $Xvar0_phy , a(cpfxcnesxyear  hospitalxdow hospitalxmonth)   resid cluster(cnes)
eststo M`h++', addscalars(Mean Mean )

qui reghdfe  m_30d  v2overlap  rw*  $XVar0 $Xvar0_phy , a(cpfxcnesxyear  hospitalxdow hospitalxmonth)   resid cluster(cnes)
eststo M`h++', addscalars(Mean Mean )

esttab M* using "${Gpath}\\table\\specialty_shares_controls.csv",   keep( *overlap) replace star(* 0.1 ** 0.05 *** 0.01) ///
se b(4) staraux brackets   ///
compress   begin(";")  delimiter(";") end(";")  scalar(Mean r2  N )

 
*************************************************************************************************
* Balance test - Physician characteristics
*************************************************************************************************
use  "${Gpath}\data\\estimation_sample.dta", clear 
keep if main_sample==1

zscore $Xvar_phy  share_age  ///
share_white share_black  share_other_race share_race_missing rais_merge overlap 

* labels 
qui do "${Gpath}\setup\\labels.do"
global NAMES ""
global NAMES_Z ""
local h=0 
estimates drop _all 
foreach i in share_age ///
share_white share_black  share_other_race share_race_missing rais_merge $Xvar_phy    {
    dis as red "`i'"
	local a=`h++'
	qui reghdfe z_`i' share_dcbo_* z_overlap  , a(cpfxcnesxyear) cluster(sp_cnes )
	estimates store M`a'
	global NAMES "${NAMES}  M`a'="${`i'}"  "
}

coefplot (M*), xline(0) keep(z_overlap) aseq swapnames ///
coeflabels( ${NAMES} , notick  labgap(2)) ciopts(recast(rcap))  eqlabels(, gap(5)) ///
graphregion(color(white) margin(zero )) ylabel(,labsize(small))  title("Panel B. Physician characteristics", pos(11) size(small)) msymbol(O) xlabel(-0.5(0.25)0.5) xlab(, nogrid) ///
xtitle("Standarized coefficients") xsize(3.2) 

graph export "${Gpath}\figures\\balance_physician.pdf", replace 

*************************************************************************************************
* Balance test - Patient characteristics
*************************************************************************************************
use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1

zscore $Xvar overlap 
* labels 
qui do "${Gpath}setup\\labels.do"
global NAMES ""
global NAMES_Z ""
local h=0 
estimates drop _all 
foreach i in $Xvar  {
    dis as red "`i'"
local a=`h++'
qui reghdfe z_`i' share_dcbo_* z_overlap  , a(cpfxcnesxyear) cluster(sp_cnes )
estimates store M`a'
global NAMES "${NAMES}  M`a'="${`i'}"  "
}

esttab M* using "${Gpath}\table\\balance_patient.csv",   keep(z_overlap) replace star(* 0.1 ** 0.05 *** 0.01) ///
se b(4) staraux brackets   compress   begin(";")  delimiter(";") end(";")  

coefplot (M*), xline(0) keep(z_overlap) aseq swapnames ///
coeflabels( ${NAMES} , notick  labgap(2)) eqlabels(, gap(5)) ///
graphregion(color(white) margin(zero )) ciopts(recast(rcap)) ///
ylabel(,labsize(small))  title("Panel A. Patient characteristics", pos(11) size(small)) msymbol(O) xlabel(-0.1(0.05)0.1) xlab(, nogrid) ///
xtitle("Standarized coefficients") xsize(3.2) ytitle("Dependent variable" " ")

graph export "${Gpath}\figures\\balance_patient.pdf", replace 


************************************************************************************************* */
* Post-study probability 
************************************************************************************************* */
use  "${Gpath}\data\\estimation_sample.dta", clear 
keep if main_sample==1 

* Baseline 
qui reghdfe  m_30d  overlap  share_dcbo_* $XVar0 $Xvar0_phy , a(cpfxcnesxyear  hospitalxdow hospitalxmonth)   resid cluster(cnes)

global df = e(df_r)
global t_crit = invttail(${df}, 0.025)  
global beta=_b[overlap]
global sebeta=_se[overlap]
global power =1-normal(${t_crit}-${beta}/${sebeta})+normal(-${t_crit}-${beta}/${sebeta})

foreach i in 0.10 0.20 0.30 {
global prior=`i'
global PSP=${power}*${prior}/(${power}*${prior}+0.05*(1-${prior}))
dis $PSP
}

************************************************************************************************* */
* Test for heterogeneous effects 
************************************************************************************************* */
use  "${Gpath}\data\\estimation_sample.dta", clear 
keep if main_sample==1 

estimates drop _all 
local h=0 

qui sum m_30d
scalar Mean=r(mean)

zscore specialty_exposure   share_shared_exp24 
* Baseline 
qui reghdfe  m_30d  share_dcbo_* $XVar0 $Xvar0_phy overlap, a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsing
eststo M`h++', addscalars(Mean Mean )

qui reghdfe  m_30d  share_dcbo_* $XVar0 $Xvar0_phy ///
z_share_shared_exp24 overlap c.(z_share_shared_exp24)#c.overlap, a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsing
eststo M`h++', addscalars(Mean Mean )

qui reghdfe  m_30d  share_dcbo_* $XVar0 $Xvar0_phy ///
z_share_shared_exp24 z_specialty_exposure overlap c.(z_share_shared_exp24 z_specialty_exposure)#c.overlap, a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsing
eststo M`h++', addscalars(Mean Mean )

esttab M* using "${Gpath}\table\\specialty_exposure.csv",   keep( z* *overlap) replace star(* 0.1 ** 0.05 *** 0.01) ///
se b(4) staraux brackets   ///
compress   begin(";")  delimiter(";") end(";")  scalar(Mean N )

************************************************************************************************* */
* Alternative samples and Same physician in the past
************************************************************************************************* */

use "${Gpath}\data\\pci_past_admission_CNS.dta", clear 

bysort SP_NAIH CNS SP_ATOPROF : keep if _n==1 
gen pastvisit=substr(SP_ATOPROF , 1, 4 )=="0301"

collapse (sum) pastvisit , by(SP_NAIH CNS) fast 
ren SP_NAIH pastN_AIH
save A.dta, replace 

use "${Gpath}data\\proceduralist_physician_level.dta", clear 
keep sp_naih CNS 
gen n_aih=sp_naih 
ren sp_naih   LN_AIH 
merge m:1 LN_AIH using "${Gpath}data\\RD_readm30.dta",  keep(1 3) gen(merge) keepusing(N_AIH)
ren N_AIH pastN_AIH
merge m:1 pastN_AIH CNS using "A.dta", gen(same_physician) keep(1 3)
erase A.dta 

collapse (max) pastvisit same_physician , by(n_aih) fast 
replace same_physician=0 if same_physician==1 
replace same_physician=1 if same_physician==3 
save A.dta, replace 

use  "${Gpath}data\\estimation_sample.dta", clear 
estimates drop _all 
local h=0 

* Baseline 
reghdfe  m_30d share_dcbo_* $XVar0 $Xvar0_phy overlap if main_sample==1, ///
a(cpfxcnesxyear hospitalxdow hospitalxmonth) cluster(sp_cnes ) keepsing
sum m_30d if e(sample)
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean )

* Exclude outliers 
qui areg  overlap share_dcbo_* if  main_sample==1 , a(cpfxcnesxyear) 
predict RES if  main_sample==1, resid 
xtile qg=RES if  main_sample==1, n(20)
gen outliers=qg==20 

reghdfe  m_30d share_dcbo_* $XVar0 $Xvar0_phy overlap if outliers==0 & main_sample==1, ///
a(cpfxcnesxyear hospitalxdow hospitalxmonth) cluster(sp_cnes ) keepsing
sum m_30d if e(sample)
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean )

* Include non-emergency patients 
qui reghdfe m_30d share_dcbo_*  $XVar0 $Xvar0_phy  overlap , a( cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes)
qui sum m_30d if e(sample)
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean)

* Exclude pandemic years 
qui reghdfe m_30d share_dcbo_*  $XVar0 $Xvar0_phy  overlap  if sp_aa<2020 & main_sample==1, a( cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes)
qui sum m_30d if e(sample)
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean)

* Exlcude proceduralist with many observations 
keep if main_sample==1 
bysort procedUnic : gen n_patients=_N 

qui reghdfe m_30d share_dcbo_*  $XVar0 $Xvar0_phy  overlap  if n_patients<2500 & main_sample==1 , ///
a( cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes)
qui sum m_30d if e(sample)
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean)

ren n_aih N_AIH 
merge m:1 N_AIH using "${Gpath}data\\RD_readm30.dta",  keep(1 3) nogen keepusing(LPROC_REA LDIAG_PRINC read*)
replace readmission30 =0 if readmission30==. 
gen subPCI=LPROC_REA=="0406030022"  | LPROC_REA=="0406030014" | LPROC_REA=="0406030030" | LPROC_REA=="0406030049"  ///
| LPROC_REA=="0406030057"  | LPROC_REA=="0406030065"  | LPROC_REA=="0406030073"
ren N_AIH  n_aih 

merge m:1 n_aih  using "A.dta", nogen  keep(3)
erase A.dta 

qui reghdfe  m_30d  share_dcbo_* $XVar0 $Xvar0_phy overlap if pastcard_intervention!=1 , a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsing
qui sum m_30d if e(sample)
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean )

esttab M* using "${Gpath}\\table\\alterative_samples.csv",  ///
keep( *overlap ) replace star(* 0.1 ** 0.05 *** 0.01) se b(4) staraux brackets   ///
compress   begin(";")  delimiter(";") end(";")  scalar(Mean N )


*******************************************************
* Inputs 
*******************************************************
use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1 

keep cpfxcnesxyear hospitalxdow hospitalxmonth  overlap sp_cnes ///
$INPUTS  val_tot lnval_tot ///
share_dcbo_* $Xvar0_phy $XVar0

zscore $INPUTS, stub(zinp_)
egen input_index=rowmean(zinp_*)

estimates drop _all
local h=0

foreach i in $INPUTS input_index val_tot  {
qui reghdfe  `i'  share_dcbo_* $Xvar0_phy $XVar0 overlap , a(cpfxcnesxyear hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsingletons
qui sum  `i' if e(sample)
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean )

}

esttab M* using "${Gpath}\\table\\Inputs.csv",   keep( overlap ) replace star(* 0.1 ** 0.05 *** 0.01) se b(4) staraux brackets   ///
compress   begin(";")  delimiter(";") end(";")  scalar(Mean N r2)

*****************************************************************************************************
* Alternative constructions of overlap 
*****************************************************************************************************
use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1
keep proc_overlap overlap  mean_overlap mainphy_overlap medi_overlap  overlap_cosine  ///
m_30d  share_dcbo_* $XVar0 $Xvar0_phy  cpfxcnesxyear  hospitalxdow hospitalxmonth sp_cnes overlap_noesp

zscore proc_overlap overlap  mean_overlap mainphy_overlap medi_overlap  overlap_cosine  overlap_noesp

estimates drop _all 
local h=0 

foreach i in  overlap  mean_overlap ///
medi_overlap mainphy_overlap  proc_overlap overlap_cosine overlap_noesp {
qui reghdfe  m_30d  share_dcbo_* $XVar0 $Xvar0_phy z_`i' , a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsing
sum m_30d if e(sample)
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean )
}

esttab M* using "${Gpath}\\table\\alterative_overlap.csv",   keep( z_*) replace star(* 0.1 ** 0.05 *** 0.01) se b(4) staraux brackets   ///
compress   begin(";")  delimiter(";") end(";")  scalar(Mean N )

*****************************************************************************************************************
* Alternative controls 
*****************************************************************************************************************
use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1

estimates drop _all
local h=1 

* no control 
qui reg m_30d overlap , cluster(sp_cnes) 
estimates store M`h++'

* proceduralist FE 
qui reghdfe m_30d share_dcbo_*  $XVar0 $Xvar0_phy overlap , a( ID_proc sp_mm  hospitalxdow hospitalxmonth) keepsing cluster(sp_cnes)
estimates store M`h++'

* hospital by day of week by month 
egen hospitalxdayofweekxmonth=group(sp_cnes day_week  sp_mm)

qui reghdfe m_30d share_dcbo_*  $XVar0 $Xvar0_phy overlap , a( cpfxcnesxyear   hospitalxdayofweekxmonth) keepsing cluster(sp_cnes)
estimates store M`h++'

* proceduralist by month and day FE
egen procxday=group(sp_cnes day_week)
egen procxmonth=group(sp_cnes sp_mm)
 
qui reghdfe m_30d share_dcbo_*  $XVar0 $Xvar0_phy overlap , a( procxday procxmonth cpfxcnesxyear) keepsing cluster(sp_cnes)
estimates store M`h++'


esttab M* using "${Gpath}\\table\\aditional_controls.csv",   keep( *overlap ) replace star(* 0.1 ** 0.05 *** 0.01) se b(4) staraux brackets   ///
compress   begin(";")  delimiter(";") end(";")  scalar(Mean N )


*********************************************************************************************************
* Controlling for diagnosis fixed effects 
*********************************************************************************************************
use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1

egen ID_cause=group(diag_princ)
estimates drop _all 
local h=1 

reghdfe  m_30d share_dcbo_* $XVar0 $Xvar0_phy overlap, a(cpfxcnesxyear hospitalxdow hospitalxmonth ID_cause) cluster(sp_cnes ) keepsing
sum m_30d if e(sample)
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean )

esttab M* using "${Gpath}\\table\\controlling_diagnosis_fixed.csv",   keep( *overlap ) replace star(* 0.1 ** 0.05 *** 0.01) se b(4) staraux brackets   ///
compress   begin(";")  delimiter(";") end(";")  scalar(Mean N )



********************************************************************************
* Nurses rate 
********************************************************************************
use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1
gen date=ym(sp_aa, sp_mm)
merge m:1 date_adm sp_cnes using  "${Gpath}data\\available_nurses_collapsed.dta", nogen keep(1 3)

bysort sp_cnes date_adm: gen npat=_N 

replace nurses_supply=1 if nurses_supply==. 

gen nure_patient=nurses_supply/npat

keep m_30d nure_patient share_dcbo_* $Xvar0_phy $XVar0 overlap cpfxcnesxyear hospitalxdow hospitalxmonth  sp_cnes

estimates drop _all 
local h=0 

qui reghdfe  nure_patient share_dcbo_* $Xvar0_phy $XVar0 overlap, a(cpfxcnesxyear hospitalxdow hospitalxmonth) cluster(sp_cnes)
qui sum  nure_patient if e(sample)
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean )

qui reghdfe m_30d  share_dcbo_* $Xvar0_phy $XVar0 overlap, a(cpfxcnesxyear hospitalxdow hospitalxmonth) cluster(sp_cnes)
qui sum  m_30d if e(sample)
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean )

qui reghdfe m_30d  nure_patient share_dcbo_* $Xvar0_phy $XVar0 overlap, a(cpfxcnesxyear hospitalxdow hospitalxmonth) cluster(sp_cnes)
qui sum  m_30d if e(sample)
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean )


esttab M* using "${Gpath}\\table\\nurses_overlap.csv",   keep(  *overlap ) replace star(* 0.1 ** 0.05 *** 0.01) ///
se b(4) staraux brackets   ///
compress   begin(";")  delimiter(";") end(";")  scalar(Mean N )



********************************************************************************
* RAIS augmented sample 
********************************************************************************
use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1

estimates drop _all 
local h=0 

qui reghdfe  m_30d  share_dcbo_* $XVar0 $Xvar0_phy overlap ///
if share_age!=. & share_white!=. , a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsing
qui sum  m_30d if e(sample)
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean )

qui reghdfe  m_30d  share_dcbo_* $XVar0 $Xvar0_phy share_age share_white  overlap, a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsing
qui sum  m_30d if e(sample)
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean )

esttab M* using "${Gpath}\\table\\RAIS_subsample.csv",   ///
keep(*overlap) replace star(* 0.1 ** 0.05 *** 0.01) se b(4) staraux brackets   ///
compress   begin(";")  delimiter(";") end(";")  scalar(Mean N )



********************************************************************************
* Total specialties 
********************************************************************************
use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1

estimates drop _all 
local h=0 

qui reghdfe  m_30d  share_dcbo_* $XVar0 $Xvar0_phy overlap , a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsing
qui sum  m_30d if e(sample)
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean )

* total specialties 
qui reghdfe  m_30d  share_dcbo_* $XVar0 $Xvar0_phy team_specialties overlap , a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsing
qui sum  m_30d if e(sample)
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean )


esttab M* using "${Gpath}\\table\\total_specilties_heart_specialties.csv",   ///
keep(team_specialties *overlap ) replace star(* 0.1 ** 0.05 *** 0.01) se b(4) staraux brackets   ///
compress   begin(";")  delimiter(";") end(";")  scalar(Mean N )


********************************************************************************
* Expertise overlap vs. shared work experience 
********************************************************************************

use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1

zscore share_shared_exp24

estimates drop _all 
local h=0 

qui reghdfe  share_shared_exp24  share_dcbo_* $Xvar0_phy $XVar0 overlap , a(cpfxcnesxyear hospitalxdow hospitalxmonth) cluster(sp_cnes)
qui sum  share_shared_exp24 if e(sample)
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean )

qui reghdfe m_30d  z_share_shared_exp24  share_dcbo_* $Xvar0_phy $XVar0 overlap , a(cpfxcnesxyear hospitalxdow hospitalxmonth) cluster(sp_cnes)
qui sum  m_30d if e(sample)
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean )


esttab M* using "${Gpath}\\table\\shared_work_threat.csv",   keep(z_share_shared_exp24  *overlap ) replace star(* 0.1 ** 0.05 *** 0.01) se b(4) staraux brackets   ///
compress   begin(";")  delimiter(";") end(";")  scalar(Mean N )


********************************************************************************
* Robust inference 
********************************************************************************


use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1

estimates drop _all 
local h=0 

*- Cluster by proceduralist 
qui reghdfe  m_30d  share_dcbo_* $XVar0 $Xvar0_phy overlap , a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(ID_proc ) keepsing
eststo M`h++'

*- Cluster by municipality 
qui reghdfe  m_30d  share_dcbo_* $XVar0 $Xvar0_phy overlap, a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(munic_mov ) keepsing
eststo M`h++'

*- Cluster by health region 
qui reghdfe  m_30d  share_dcbo_* $XVar0 $Xvar0_phy overlap, a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(co_regsaud_mov) keepsing
eststo M`h++'


*two-way cluster 
reghdfe m_30d share_dcbo_*  $XVar0 $Xvar0_phy overlap, a(  cpfxcnesxyear hospitalxdow hospitalxmonth) cluster(date sp_cnes) keepsingletons
eststo M`h++'

esttab M* using "${Gpath}\\table\\robust_inference.csv",   keep( *overlap) ///
replace star(* 0.1 ** 0.05 *** 0.01) se b(4) staraux brackets   ///
compress   begin(";")  delimiter(";") end(";")   scalar( p_wild se_young N r2)

*****************************************************************************************************************
* LASSO
*****************************************************************************************************************
use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1

estimates drop _all 
local h=0 

* Baseline 
qui reghdfe  m_30d  share_dcbo_* $XVar0 $Xvar0_phy overlap, a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsing
eststo M`h++'

*-- patient 
qui rlassologit m_30d $XVar0 
local VAR =e(selected)

global LASSO_1 ""
foreach i in `VAR' {    
global LASSO_1 "${LASSO_1} `i'"	
}

qui reghdfe  m_30d  share_dcbo_* ${LASSO_1} overlap, a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsing
eststo M`h++'


*-- patient +interactions 
qui rlassologit m_30d $XVar0  ${INTPX}
local VAR =e(selected)

global LASSO_2 ""
foreach i in `VAR' {    
global LASSO_2 "${LASSO_2} `i'"	
}

qui reghdfe  m_30d  share_dcbo_* ${LASSO_2} overlap, a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsing
eststo M`h++'

*-- physician 
qui rlassologit m_30d $Xvar0_phy 
local VAR =e(selected)

global LASSO_3 ""
foreach i in `VAR' {    
global LASSO_3 "${LASSO_3} `i'"	
}

qui reghdfe  m_30d  share_dcbo_* ${LASSO_3} overlap, a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsing
eststo M`h++'

*-- physician + interactions 
qui rlassologit m_30d $Xvar0_phy ${INTDX} 
local VAR =e(selected)

global LASSO_4 ""
foreach i in `VAR' {    
global LASSO_4 "${LASSO_4} `i'"	
}

qui reghdfe  m_30d  share_dcbo_* ${LASSO_4} overlap, a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsing
eststo M`h++'


* Patient + physician 
qui rlassologit m_30d $Xvar0_phy ${INTDX} $XVar0  ${INTPX}
local VAR =e(selected)

global LASSO_5 ""
foreach i in `VAR' {    
global LASSO_5 "${LASSO_5} `i'"	
}

qui reghdfe  m_30d  share_dcbo_* ${LASSO_5} overlap, a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsing
eststo M`h++'

esttab M* using "${Gpath}\\table\\LASSO.csv",   keep( overlap ) replace star(* 0.1 ** 0.05 *** 0.01) se b(4) staraux brackets   ///
compress   begin(";")  delimiter(";") end(";")  scalar(Mean N )

**************************************************************************************************
* Emily Oster 
**************************************************************************************************
use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1

gen B_0=. 
gen R2_0=. 
gen R2=. 
gen B_1=. 
gen Delta=. 
gen model=_n 
replace model=. if model>2 

* Basleine 
qui reghdfe m_30d share_dcbo_* overlap, a( cpfxcnesxyear)
replace B_0=_b[overlap] if model!=. 
replace R2_0=e(r2) if model!=. 

* add patient and physician characteristics
qui reghdfe m_30d share_dcbo_* past_read $XVar0 $Xvar0_phy overlap, a( cpfxcnesxyear )
replace B_1=_b[overlap] if model==1 
replace R2=e(r2) if model==1 

replace Delta=((B_1)/(B_0-B_1))*( (R2-R2_0)/(0.3*R2) )

mkmat B_1 B_0 R2 R2_0 Delta if model!=. , mat(A)
esttab matrix(A) using "${Gpath}\\table\\Oster_test.csv", replace 



*************************************************************************************************
* Exclude specialty configurartions 
*************************************************************************************************

use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1

keep coverlap* hospitalxdow hospitalxmonth sp_cnes m_30d  share_dcbo_* $XVar0 $Xvar0_phy cpfxcnesxyear 
gen beta=. 
gen sebeta=.
gen id=_n 
forvalues i=1/50 {
dis as red "reg `i'"
qui reghdfe  m_30d  share_dcbo_* $XVar0 $Xvar0_phy coverlap`i', a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsing
replace beta=_b[coverlap`i'] if id==`i'
replace sebeta=_se[coverlap`i'] if id==`i' 
}

keep beta sebeta id 
drop if beta==. 
saveold "${Gpath}data\\leave_one_out_results.dta", replace 

use "${Gpath}data\\leave_one_out_results.dta", clear 
erase "${Gpath}data\\leave_one_out_results.dta"
gen lb=beta-1.95*sebeta
gen hb=beta+1.95*sebeta

twoway (rcap lb hb id) (scatter beta id), graphregion(color(white) margin(0)) ylabel(-0.08(0.04)0) ///
xlabel(1(2)50) legend(off) xtitle(" " "Proceduralist-physician specialty configuration excluded") ///
ytitle(" ") xsc(range(1 52))

graph export "${Gpath}\figures\\leave_one_out_results.pdf", replace 

************************************************************************************************* */
* Selection into procedure 
************************************************************************************************* */

use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1

merge 1:1 n_aih using "${Gpath}data\\pci_propensity.dta", keep(1 3) nogen 

estimates drop _all 
local h=0 
* Baseline specification 
qui areg predict_pci  share_dcbo_* overlap , a(cpfxcnesxyear) cluster(sp_cnes)
qui sum  predict_pci if e(sample)
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean )

* hospital x (month- and day-of-week) FE 
qui reghdfe  predict_pci  share_dcbo_* overlap, a(cpfxcnesxyear hospitalxdow hospitalxmonth) cluster(sp_cnes)
qui sum  predict_pci if e(sample)
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean )

* patient characteristics
qui areg predict_pci  share_dcbo_* overlap $XVar0 , a(cpfxcnesxyear) cluster(sp_cnes)
qui sum  predict_pci if e(sample)
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean )

* Physician characteristics
qui areg predict_pci  share_dcbo_* overlap $Xvar0_phy , a(cpfxcnesxyear) cluster(sp_cnes)
qui sum  predict_pci if e(sample)
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean )

* All 
qui reghdfe  predict_pci  share_dcbo_* overlap $Xvar0_phy $XVar0 , a(cpfxcnesxyear hospitalxdow hospitalxmonth) cluster(sp_cnes)
qui sum  predict_pci if e(sample)
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean )


esttab M* using "${Gpath}table\\predicted_PCI_actual_overlap.csv",   keep( *overlap ) replace star(* 0.1 ** 0.05 *** 0.01) se b(6) staraux brackets ///
compress   begin(";")  delimiter(";") end(";")  scalar(Mean N )


estimates drop _all 
local h=0 

qui reghdfe  m_30d share_dcbo_* $XVar0 $Xvar0_phy overlap  , a(cpfxcnesxyear hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsin
qui sum m_30d if e(sample) 
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean)

qui reghdfe  m_30d share_dcbo_* $XVar0 $Xvar0_phy overlap  predict_pci , a(cpfxcnesxyear hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsin
qui sum m_30d if e(sample) 
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean)

qui reghdfe  m_30d share_dcbo_* $XVar0 $Xvar0_phy overlap predict_pci c.predict_pci#c.predict_pci  , a(cpfxcnesxyear hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsin
qui sum m_30d if e(sample) 
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean)

qui reghdfe  m_30d share_dcbo_* $XVar0 $Xvar0_phy overlap  predict_pci c.predict_pci#c.predict_pci ///
c.predict_pci#c.predict_pci#c.predict_pci , a(cpfxcnesxyear hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsin
qui sum m_30d if e(sample) 
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean)

qui reghdfe  m_30d share_dcbo_* $XVar0 $Xvar0_phy overlap predict_pci c.predict_pci#c.predict_pci ///
c.predict_pci#c.predict_pci#c.predict_pci c.predict_pci#c.predict_pci#c.predict_pci#c.predict_pci , a(cpfxcnesxyear hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsin
qui sum m_30d if e(sample) 
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean)

qui reghdfe  m_30d share_dcbo_* $XVar0 $Xvar0_phy overlap  mills  , a(cpfxcnesxyear hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsin
qui sum m_30d if e(sample) 
scalar Mean=r(mean)
eststo M`h++', addscalars(Mean Mean)



esttab M* using "${Gpath}\\table\\selection_into_PCI.csv",   keep( *overlap ) replace star(* 0.1 ** 0.05 *** 0.01) se b(4) staraux brackets   ///
compress   begin(";")  delimiter(";") end(";")  scalar(Mean N )


********************************************************************************
* Summary statistics
********************************************************************************
use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1

cap erase "${Gpath}\\table\\summary_statistics_main.csv"
gen age_65_79=idade>=65 & idade<=79
gen age_65_less=idade<65

gen heart_attack=substr(diag_princ , 1, 3)=="I21" | substr(diag_princ , 1, 3)=="I22" | substr(diag_princ , 1, 3)=="I23"
gen angina_pectoris=substr(diag_princ , 1, 3)=="I20"
gen other_acute=substr(diag_princ , 1, 3)=="I24"
gen chronic_heart_isch=substr(diag_princ , 1, 3)=="I25"
gen other_disease=1-chronic_heart_isch-other_acute-angina_pectoris-heart_attack


qui sum overlap 
mat A=(r(mean), r(sd))

global NAMES "overlap"
foreach i in m_30d idade age_over80 age_65_79 age_65_less male race_1 race_2 race_other ///
heart_attack angina_pectoris other_acute chronic_heart_isch other_disease  lstay val_tot  {
qui sum `i'
mat A=A\(r(mean), r(sd))
global NAMES "${NAMES} `i'"
}

global N=_N 
mat A=A\(., ${N})

qui unique sp_cnes
mat A=A\(., r(unique) )

qui unique procedUnic
mat A=A\(., r(unique) )

ren n_aih sp_naih
merge 1:m sp_naih using "${Gpath}data\\proceduralist_physician_level.dta",  keepusing(proceduralist CNS) keep(3) nogen  
drop if proceduralist==1

qui unique CNS 
mat A=A\(., r(unique) )

matrix rownames A=${NAMES} n_patients n_hospitals n_proceduralist n_physicians

esttab matrix(A) using "${Gpath}\\table\\summary_statistics_main.csv",   replace ///
compress   begin(";")  delimiter(";") end(";")  

********************************************************************************
* Variation in expertise
********************************************************************************
use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1

areg overlap share_dcbo_*, a(cpfxcnesxyear)
predict resx0, resid


hist resx0 , frac  xtitle("Expertise overlap" "(residualized)") ///
title("Panel A. Histogram of the residualized expertise overlap", pos(11) size(medium))  graphregion(color(white) margin(zero ))
graph export "${Gpath}\figures\\histogram_residual.pdf", replace 

qui reg overlap share_dcbo_*
predict resx, resid
bysort cpfxcnesxyear: egen mean=mean(resx)
gen devmedian=resx-mean
gen devmean=resx-mean

collapse (iqr) iqr=resx (mean) devmean (median) devmedian (min) min=resx (max) max=resx (sd) sd=resx, by(cpfxcnesxyear )
gen maxmin=max-min
hist iqr , frac

hist sd ,   frac  xtitle("Standard deviation in residualized expertise overlap") ///
title("Panel B. Standard deviation within proceduralist", pos(11) size(medium))  graphregion(color(white) margin(zero ))
graph export "${Gpath}\figures\\histogram_residual_SD.pdf", replace 

hist maxmin ,  frac  xtitle("Max-min difference in residualized expertise overlap") ///
title("Panel C. Max-min difference  within proceduralist", pos(11) size(medium))  graphregion(color(white) margin(zero ))
graph export "${Gpath}\figures\\histogram_residual_MINMAX.pdf", replace 

hist iqr ,   frac  xtitle("Interquartile range in residualized expertise overlap") ///
title("Panel D. Interquartile range within proceduralist", pos(11) size(medium))  graphregion(color(white) margin(zero ))
graph export "${Gpath}\figures\\histogram_residual_iqr.pdf", replace 

********************************************************************************
* Spending histogram 
********************************************************************************
use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1

hist val_tot  , /*fcolor(gs9) lcolor(black)*/   frac  xtitle("") ///
title(" ", pos(11) size(medium))  graphregion(color(white) margin(zero )) ///
xlabel(0 "0" 20000 "20000" 40000 "40000" 60000 "60000" 80000 "80000" 100000 "100000" )

graph export "${Gpath}\figures\\spending_hospitals.pdf", replace 


********************************************************************************************************************
* Number of patient per Hospitals 
********************************************************************************************************************

use  "${Gpath}data\\estimation_sample.dta", clear 

keep if main_sample==1 

keep sp_cnes
bysort sp_cnes: gen n_patients=_N 
duplicates drop n_patients, force 
sort n_patients

sum n_patients, detail 

gen share=n_patients/r(sum)
gen cum_shre=sum(share)

local N=_N 
local a=round(0.1*`N')
dis cum_shre[`N']-cum_shre[`N'-`a'-1]

hist n_patients ,    frac  xtitle("") ///
title("Panel A. Number of patients per hospital", pos(11) size(medium))  graphregion(color(white) margin(zero ))

graph export "${Gpath}\figures\\patients_hospitals.pdf", replace 


********************************************************************************************************************
* Number of patient per proceduralist
********************************************************************************************************************

use  "${Gpath}data\\estimation_sample.dta", clear 

keep if main_sample==1 

keep procedUnic 
bysort procedUnic : gen n_patients=_N 
duplicates drop n_patients, force 
sort n_patients

sum n_patients, detail 

gen share=n_patients/r(sum)
gen cum_shre=sum(share)

local N=_N 
local a=round(0.1*`N')
dis cum_shre[`N']-cum_shre[`N'-`a'-1]


hist n_patients ,    frac  xtitle("") ///
title("Panel B. Number of patients per proceduralist", pos(11) size(medium))  graphregion(color(white) margin(zero ))

graph export "${Gpath}\figures\\patients_proceduralist.pdf", replace 




*******************************************************
* Heterogeneity by shared work experience 
*******************************************************

use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1 

qui sum share_shared_exp24, detail  
local a25=r(p25)
local a50=r(p50)
local a75=r(p75)

sum m_30d if share_shared_exp24<`a25' & share_shared_exp24!=. 
local m25=r(mean)
local sd25=r(sd)

sum m_30d if share_shared_exp24>`a75' & share_shared_exp24!=. 
local m75=r(mean)
local sd75=r(sd)

sum m_30d if share_shared_exp24>`a25' & share_shared_exp24<`a75' & share_shared_exp24!=. 
local m50=r(mean)
local sd50=r(sd)

gen y0=m_30d
replace y0=(m_30d/`m25')  if share_shared_exp24<`a25' & share_shared_exp24!=. 
replace y0=(m_30d/`m50')  if share_shared_exp24>`a25' & share_shared_exp24<`a75' & share_shared_exp24!=. 
replace y0=(m_30d/`m75')  if share_shared_exp24>`a75' & share_shared_exp24!=. 

* Baseline 
reghdfe m_30d  share_dcbo_*  $Xvar0_phy $XVar0  c.share_shared_exp24#c.share_shared_exp24  overlap ///
share_shared_exp24 c.overlap#c.share_shared_exp24 ///
c.overlap#c.share_shared_exp24#c.share_shared_exp24 , a(cpfxcnesxyear hospitalxdow hospitalxmonth) cluster(sp_cnes)

gen percentile=_n 
gen beta=. 
gen sebeta=. 

foreach i in 25 50 75 {
local c1=`a`i''
local c2=`a`i''^2
	
lincom overlap +c.overlap#c.share_shared_exp24*`c1' +c.overlap#c.share_shared_exp24#c.share_shared_exp24*`c2'
replace beta=r(estimate) if percentile==`i'
replace sebeta=r(se) if percentile==`i'
}

gen low=beta-1.96*sebeta
gen high=beta+1.96*sebeta

preserve 
drop if low==. 
twoway (bar beta percentile if percentile==25,  lcolor(navy)  barwidth(20) ) ///
(bar beta percentile if percentile==50,   lcolor(maroon) barwidth(20) ) ///
(bar beta percentile if percentile==75, lcolor(green)   barwidth(20) ) ///
(rcap low high percentile, lcolor(black) ) , ylabel(-0.06(0.02)0.02) graphregion(color(white) margin(zero)) legend(off) ///
xlabel(25 "25th percentile" 50 "50t percentile" 75 "75th percentile") ///
xtitle("Percentile of the shared work experience distribution", color(blue)) ytitle("Effect of expertise overlap", color(blue)) ///
title("Panel B. Low versus High Shared Work Experience", pos(11))

graph export "${Gpath}\figures\\heterogeneity_experience_together.pdf", replace 
restore 

foreach i in low high beta sebeta percentile {
    drop `i'
}

* Robustness 
reghdfe y0 share_dcbo_*  $Xvar0_phy $XVar0  c.share_shared_exp24#c.share_shared_exp24  overlap ///
share_shared_exp24 c.overlap#c.share_shared_exp24 ///
c.overlap#c.share_shared_exp24#c.share_shared_exp24 , a(cpfxcnesxyear hospitalxdow hospitalxmonth) cluster(sp_cnes)

gen percentile=_n 
gen beta=. 
gen sebeta=. 

foreach i in 25 50 75 {
local c1=`a`i''
local c2=`a`i''^2
	
lincom overlap +c.overlap#c.share_shared_exp24*`c1' +c.overlap#c.share_shared_exp24#c.share_shared_exp24*`c2'
replace beta=(r(estimate)) if percentile==`i'
replace sebeta=(r(se)) if percentile==`i'
}

gen low=beta-1.96*sebeta
gen high=beta+1.96*sebeta

preserve 
drop if low==. 
twoway (bar beta percentile if percentile==25,  lcolor(navy)  barwidth(20) ) ///
(bar beta percentile if percentile==50,   lcolor(maroon) barwidth(20) ) ///
(bar beta percentile if percentile==75, lcolor(green)   barwidth(20) ) ///
(rcap low high percentile, lcolor(black) ) , ylabel(-1(0.5)0.5) graphregion(color(white) margin(zero)) legend(off) ///
xlabel(25 "25th percentile" 50 "50t percentile" 75 "75th percentile") ///
xtitle("Percentile of the shared work experience distribution", color(blue)) ytitle("Effect of expertise overlap", color(blue)) ///
title("Panel B. Low versus High Shared Work Experience", pos(11))

graph export "${Gpath}\figures\\heterogeneity_experience_together_robustness.pdf", replace 
restore 

*******************************************************
* Heterogeneity by mortality risk 
*******************************************************
use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1

* predicted mortality  
qui do "${Gpath}\setup\\interactions.do"
cap drop dx*

gen cause_3d=substr(diag_princ, 1, 3)
egen cause_code=group(cause_3d)

probit m_30d $XVar0 ${INTPX}   i.cause_code
predict pr_y1, pr 
drop px*

gen y0=m_30d

sum pr_y1, detail  
local a25=r(p25)
local a50=r(p50)
local a75=r(p75)

sum m_30d if pr_y1<`a25' & pr_y1!=. 
local m25=r(mean)
local sd25=r(sd)

sum m_30d if pr_y1>`a75' & pr_y1!=. 
local m75=r(mean)
local sd75=r(sd)

sum m_30d if pr_y1>`a25' & pr_y1<`a75' & pr_y1!=. 
local m50=r(mean)
local sd50=r(sd)



replace y0=(m_30d/`m25') if pr_y1<`a25' & pr_y1!=. 
replace y0=(m_30d/`m50') if pr_y1>`a25' & pr_y1<`a75' & pr_y1!=. 
replace y0=(m_30d/`m75') if pr_y1>`a75' & pr_y1!=. 

* Baseline 
reghdfe m_30d  share_dcbo_*  $Xvar0_phy $XVar0 pr_y1 c.pr_y1#c.pr_y1  overlap ///
c.overlap#c.pr_y1 c.overlap#c.pr_y1#c.pr_y1  , a(cpfxcnesxyear hospitalxdow hospitalxmonth) cluster(sp_cnes)

test c.overlap#c.pr_y1 c.overlap#c.pr_y1#c.pr_y1
global F=r(F)
global F=substr("${F}", 1, 5)

cap drop percentile 
cap drop beta 
cap drop sebeta 
gen percentile=_n 
gen beta=. 
gen sebeta=. 

foreach i in 25 50 75 {

local c1=`a`i''	
local c2=`a`i''^2
lincom overlap+c.overlap#c.pr_y1*`c1' +c.overlap#c.pr_y1#c.pr_y1*`c2' 
replace beta=r(estimate) if percentile==`i'
replace sebeta=r(se) if percentile==`i'
}


cap drop low 
cap drop high 
gen low=beta-1.96*sebeta
gen high=beta+1.96*sebeta

preserve 
drop if low==. 
twoway (bar beta percentile if percentile==25,  lcolor(navy)  barwidth(20) ) ///
(bar beta percentile if percentile==50,   lcolor(maroon) barwidth(20) ) ///
(bar beta percentile if percentile==75, lcolor(green)   barwidth(20) ) ///
(rcap low high percentile, lcolor(black) ) , ylabel(-0.08(0.02)0.02) graphregion(color(white) margin(zero)) legend(off) ///
xlabel(25 "25th percentile" 50 "50th percentile" 75 "75th percentile") ///
xtitle("Percentile of the predicted mortality distribution", color(blue)) ytitle("Effect of expertise overlap", color(blue)) ///
title("Panel A. Low versus High Predicted Mortality Risk", pos(11))

graph export "${Gpath}\figures\\heterogeneity_mortality_risk.pdf", replace 
restore 

foreach i in low high beta sebeta percentile {
    drop `i'
}

* Robustness
reghdfe y0 share_dcbo_*  $Xvar0_phy $XVar0 pr_y1 c.pr_y1#c.pr_y1  overlap ///
c.overlap#c.pr_y1 c.overlap#c.pr_y1#c.pr_y1  , a(cpfxcnesxyear hospitalxdow hospitalxmonth) cluster(sp_cnes)

test c.overlap#c.pr_y1 c.overlap#c.pr_y1#c.pr_y1
global F=r(F)
global F=substr("${F}", 1, 5)

cap drop percentile 
cap drop beta 
cap drop sebeta 
gen percentile=_n 
gen beta=. 
gen sebeta=. 

foreach i in 25 50 75 {

local c1=`a`i''	
local c2=`a`i''^2
lincom overlap+c.overlap#c.pr_y1*`c1' +c.overlap#c.pr_y1#c.pr_y1*`c2' 
replace beta= (r(estimate)) if percentile==`i'
replace sebeta=(r(se)) if percentile==`i'

}


cap drop low 
cap drop high 
gen low=beta-1.96*sebeta
gen high=beta+1.96*sebeta

preserve 
drop if low==. 
twoway (bar beta percentile if percentile==25,  lcolor(navy)  barwidth(20) ) ///
(bar beta percentile if percentile==50,   lcolor(maroon) barwidth(20) ) ///
(bar beta percentile if percentile==75, lcolor(green)   barwidth(20) ) ///
(rcap low high percentile, lcolor(black) ) , ylabel(-1(0.5)0.5, nogrid) graphregion(color(white) margin(zero)) legend(off) ///
xlabel(25 "25th percentile" 50 "50th percentile" 75 "75th percentile") ///
xtitle("Percentile of the predicted mortality distribution", color(blue)) ytitle("Effect of expertise overlap", color(blue)) ///
title("Panel A. Low versus High Predicted Mortality Risk", pos(11))

graph export "${Gpath}\figures\\heterogeneity_mortality_risk_robustness.pdf", replace 
restore 

*************************************************************************************************
* Effects of Past, Actual and Future Expertise Overlap
*************************************************************************************************
*-- process 
use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1
keep sp_aa sp_mm  overlap procedUnic

collapse (mean) Exp_wmexp=overlap, by(sp_aa sp_mm procedUnic)

destring sp_aa sp_mm, replace 
gen date =ym(sp_aa, sp_mm)
egen id=group(procedUnic)
xtset id date , monthly 

forvalues i=1/5 {
gen L`i'Exp_wmexp=L`i'.Exp_wmexp
gen F`i'Exp_wmexp=F`i'.Exp_wmexp
}

keep sp_aa sp_mm procedUnic F* L* 

* test for serial correlation 
areg L1Exp_wmexp i.sp_aa  i.sp_mm L2Exp_wmexp, a(procedUnic )

tempfile leads_lags_exp
saveold   "`leads_lags_exp'.dta", replace 

*-- Estimation 
use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1

merge m:1 sp_aa sp_mm  procedUnic using "`leads_lags_exp'.dta", nogen keep(1 3)

gen Beta=. 
gen seBeta=. 

gen case_time=_n 
replace case_time=case_time-6
replace case_time=. if case_time>5 

forvalues i=1(1)4 {
dis as red "lag `i'"
qui reghdfe  m_30d $XVar0 $Xvar_phy  share_dcbo_*  L`i'Exp_wmexp  , a(cpfxcnesxyear hospitalxdow hospitalxmonth) cluster(sp_cnes)
replace Beta=_b[L`i'Exp_wmexp] if case_time==-`i'
replace seBeta=_se[L`i'Exp_wmexp] if case_time==-`i'

dis as red "Lead `i'"

qui reghdfe   m_30d $XVar0 $Xvar_phy  share_dcbo_*  F`i'Exp_wmexp , a(cpfxcnesxyear hospitalxdow hospitalxmonth) cluster(sp_cnes)
replace Beta=_b[F`i'Exp_wmexp] if case_time==`i'
replace seBeta=_se[F`i'Exp_wmexp] if case_time==`i'
}

qui reghdfe m_30d share_dcbo_* $Xvar0_phy $XVar0 overlap, a(cpfxcnesxyear hospitalxdow hospitalxmonth) cluster(sp_cnes)
replace Beta=_b[overlap] if case_time==0
replace seBeta=_se[overlap] if case_time==0

gen Lbeta=Beta-1.95*seBeta
gen Hbeta=Beta+1.95*seBeta
 
gen zeroline=0

cap drop y1 y2 
gen y1=0.04
gen y2=-0.07


twoway (rarea y1 y2 case_time if case_time<=-1 & case_time>=-4 , color(gs15) ) (rarea y1 y2 case_time if case_time>=1 & case_time<=4, color(gs15)) ///
(line zeroline case_time if abs(case_time)<=4, lp(dash) lcolor(black))  ///
(scatter Lbeta Hbeta case_time  if abs(case_time)<=4  , lw(vthin vthin) connect(l l) lcolor(black black) lp(dash dash) msymbol(none none) )   ///
(scatter Beta case_time if abs(case_time)<=4, connect(l) lcolor(black) mcolor(black) lw(medthick) msymbol(O)) ////
, graphregion(color(white)) ///
text(-0.04 -2.5 "Placebo:" "Past expertise overlap") ///
text(-0.04 2.5 "Placebo:" "Future expertise overlap") ///
text(0.02 0 "Actual" "expertise overlap" ) ///
ylabel(-0.06(0.02)0.04, nogrid) xlabel(-4(1)4) legend(off) xtitle("Months relative to actual case") xsize(6)

graph export "${Gpath}figures\\lags_leads_effects.pdf", replace 


*************************************************************************************************
* Baseline - Visual evidence 
*************************************************************************************************
use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1

areg m_30d share_dcbo_* , a(cpfxcnesxyear )
predict res_y, resid 
qui sum m_30d
replace res_y=r(mean)+res_y

areg overlap share_dcbo_*  , a(cpfxcnesxyear )
predict res_x, resid 

qui sum res_x , detail 
gen above_median=res_x>=r(p50)
gen below_median=1-above_median
gen treatment=above_median

qui areg m_30d share_dcbo_* overlap, a(cpfxcnesxyear)  cluster(sp_cnes)

global b=string(round(_b[overlap], 0.0001),"%9.3f")
global se=string(round(_se[overlap], 0.0001),"%9.3f")

binscatter res_y res_x , n(20)   ytitle("") ///
ytitle("Actual 30-day mortality") xtitle("Expertise overlap" "(residualized)") title(" ", size(medium) pos(11) )  ///
text( 0.025 0.2 "coeff.= ${b}" "   s.e.= (${se})") ///
graphregion(color(white) margin(zero )) ylabel(0(0.02)0.1) xlabel(-0.3(0.1)0.3)
graph export "${Gpath}\figures\\actual_mortality_overlap.pdf", replace 


*************************************************************************************************
* Identification startegy 
*************************************************************************************************
clear 
set obs 4 
gen group=_n 
gen z=1 if group==1 
replace z=0.01 if group==2 
replace z=0.5 if group>=3 

twoway (bar z group, barw(0.65) color(none)  lw(medium) ) ///
(bar z group if group==1 | group==3, barw(0.65) color(navy) fi(120) lw(medium) ) ///
(bar z group if group==2 | group==4, barw(0.65) color(navy) fi(50) lw(medium) ) ///
(pci 1 2 0.043 2  1 2 1 1.9 0.043 2 0.043 1.9, lc(black)) ///
, ylabel(0(0.5)1, nogrid) ///
xline(2.5, lcolor(black) lp(dash)) xlabel(1 `"  "Cardiology"    "'  2  "Internal Medicine"  ///
3 `"   "Cardiology"    "'  4  "Internal Medicine" ) xtitle(" ""Physician's specialty") ///
ytitle("Expertise overlap") ///
text(1.35 1.5 "Proceduralist A") text(1.35 3.5 "Proceduralist B") ///
text(1.25 1.5   "Cardiology"  , size(small)) ///
text(1.2 3.5  "Cardiology" "+" "Internal Medicine"  , size(small)) ///
text( 0.5 2.1 "1") ///
title("", pos(11) ) ///
legend(off) graphregion(color(white) margin(0)) plotregion(ilcolor(white)) ysc(range(0 1.4)) xsc(range(1 4.4))

graph export "${Gpath}\figures\\identification_A.pdf", replace 


*************************************************************************************************
* Selected combinations of specialties -- Predicted mortality 
*************************************************************************************************

use "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample ==1 

qui logit m_30d $XVar0 
predict pr_y1, pr 

keep n_aih m_30d pr_y1
tempfile A 
save `A'.dta, replace 

use "${Gpath}data\\proceduralist_physician_level.dta",  clear 

keep code1 code2 cod_olap  sp_naih tot_vi overlap
ren sp_naih n_aih 
merge m:1 n_aih  using "`A'.dta", nogen keepusing(n_aih pr_y1) keep(3)
gen cases=1 
set seed 123 
gen random=runiform()
gsort n_aih -tot_vi random
by n_aih: keep if _n==1 
collapse (mean) overlap m_30d =pr_y1 (sum) cases  , by(code1 code2 ) 


keep if (code1=="8 17" & code2=="8 17" ) ///
| (code1=="8 17" & code2=="8" ) ///
| (code1=="17 43" & code2=="8" ) ///
| (code1=="8 25" & code2=="8 17" ) ///
| (code1=="8 9 12 17" & code2=="8 9 17" ) ///
| (code1=="8 17 19"  & code2=="8 17 25") ///
| (code1=="8 17 25"  & code2=="8")

gsort -overlap -m_30d
gen g=8-_n 

gen x =real(substr(string(m_30d), 1, 4))
drop m_30d
ren x m_30d

twoway (scatter m_30d g, color(none) xaxis(1 2)) ///
(bar m_30d g, color(navy)) ///
(pcarrowi 0.097 2.5  0.097 1.3, color(black)) ///
(pcarrowi 0.097 5.5 0.097 6.7  , color(black)),  ylabel(0.02(0.03)0.11, nogrid) ///
xlabel(1 `""Internal Med." "+"   "Nephrology" "' /// 
2 `""Cardiology." "+" "Critical Care Med." "' ///
3 `""Internal Med."  "+" "Cardiology" "+" "Critical Care Med." "' ///
4 `""Internal Med" "+" "Cardiology" "+" "Fam. Med.""' ///
5 `""Internal Med." "+"   "Cardiology" "' ///
6 `" "Internal Med." "+" "Cardiology" "+" "General Surg." "+" "Cardiovasc. Surg." "' ///
7 `""Internal Med." "+"  "Cardiology"   "' ///
, axis(1) labsize(small)) ///
xsize(6.7) xtitle("{bf:Proceduralist}" " ", axis(2)) xtitle(" " "{bf:Physician} ", axis(1)) legend(off) ///
xlabel(1  "Cardiology" /// 
2 `""Internal Med." "+" "Cardiology""' ///
3 "Cardiology" ///
4  `""Internal Med." "+" "Cardiology" "+"  "Critical Care Med.""' ///
5 `"Cardiology "' ///
6  `""Internal Med." "+" "Cardiology" "+" "Cardiovasc. Surg."  "' ///
7 `""Internal Med." "+"  "Cardiology" "' ///
, axis(2) labsize(small)) ytitle("30-day mortality") graphregion(color(white)) ///
text(0.105 6 "higher overlap") ///
text(0.105 2 "lower overlap")

graph export "${Gpath}\figures\\proceduralist_physician_comparisons_combined_pr.pdf", replace 


*************************************************************************************************
* Selected combinations of specialties 
*************************************************************************************************

use "${Gpath}data\\proceduralist_physician_level.dta",  clear 

keep code1 code2 cod_olap  sp_naih tot_vi overlap
ren sp_naih n_aih 
merge m:1 n_aih  using "${Gpath}data\\estimation_sample.dta", nogen keepusing(n_aih m_30d main_sample) keep(3)
keep if main_sample==1 
gen cases=1 
set seed 123 
gen random=runiform()
gsort n_aih -tot_vi random
by n_aih: keep if _n==1 
collapse (mean) overlap m_30d (sum) cases  , by(code1 code2 ) 


keep if (code1=="8 17" & code2=="8 17" ) ///
| (code1=="8 17" & code2=="8" ) ///
| (code1=="17 43" & code2=="8" ) ///
| (code1=="8 25" & code2=="8 17" ) ///
| (code1=="8 9 12 17" & code2=="8 9 17" ) ///
| (code1=="8 17 19"  & code2=="8 17 25") ///
| (code1=="8 17 25"  & code2=="8")

gsort -overlap -m_30d
gen g=8-_n 

gen x =real(substr(string(m_30d), 1, 4))
drop m_30d
ren x m_30d

twoway (scatter m_30d g, color(none) xaxis(1 2)) ///
(bar m_30d g, color(navy)) ///
(pcarrowi 0.097 2.5  0.097 1.3, color(black)) ///
(pcarrowi 0.097 5.5 0.097 6.7  , color(black)),  ylabel(0.02(0.03)0.11, nogrid) ///
xlabel(1 `""Internal Med." "+"   "Nephrology" "' /// 
2 `""Cardiology." "+" "Critical Care Med." "' ///
3 `""Internal Med."  "+" "Cardiology" "+" "Critical Care Med." "' ///
4 `""Internal Med" "+" "Cardiology" "+" "Fam. Med.""' ///
5 `""Internal Med." "+"   "Cardiology" "' ///
6 `" "Internal Med." "+" "Cardiology" "+" "General Surg." "+" "Cardiovasc. Surg." "' ///
7 `""Internal Med." "+"  "Cardiology"   "' ///
, axis(1) labsize(small)) ///
xsize(6.7) xtitle("{bf:Proceduralist}" " ", axis(2)) xtitle(" " "{bf:Physician} ", axis(1)) legend(off) ///
xlabel(1  "Cardiology" /// 
2 `""Internal Med." "+" "Cardiology""' ///
3 "Cardiology" ///
4  `""Internal Med." "+" "Cardiology" "+"  "Critical Care Med.""' ///
5 `"Cardiology "' ///
6  `""Internal Med." "+" "Cardiology" "+" "Cardiovasc. Surg."  "' ///
7 `""Internal Med." "+"  "Cardiology" "' ///
, axis(2) labsize(small)) ytitle("30-day mortality") graphregion(color(white)) ///
text(0.105 6 "higher overlap") ///
text(0.105 2 "lower overlap")

graph export "${Gpath}\figures\\proceduralist_physician_comparisons_combined.pdf", replace 


 
*************************************************************************************************
* Raw correlation between expertise overlap and mortality 
*************************************************************************************************
use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1

binscatter m_30d overlap , xtitle("Expertise overlap") ///
ytitle("30-day mortality") graphregion(color(white) margin(0)) n(20) ylabel(0(0.02)0.1)

graph export "${Gpath}\figures\\raw_correlation_expertise_m30day.pdf", replace 


*****************************************************************************************************
* Demographic concordance and Expertise Overlap
*****************************************************************************************************
use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1
zscore proc_overlap overlap  mean_overlap mainphy_overlap medi_overlap  overlap_cosine sim_age same_race same_gender 

gen z_demo_similarity =z_sim_age +z_same_race +z_same_gender 

binscatter same_gender overlap, ylabel(0(0.2)1) ytitle("Same gender") xtitle("Expertise overlap") n(20) graphregion(color(white) margin(0)) 
graph export "${Gpath}\figures\\same_gender_overlap.pdf", replace 

binscatter same_race overlap, ylabel(0(0.2)1) ytitle("Same race") xtitle("Expertise overlap") n(20) graphregion(color(white) margin(0)) 
graph export "${Gpath}\figures\\same_race_overlap.pdf", replace 

gen age_dist=-sim_age 
binscatter age_dist overlap, ylabel(0(10)30) ytitle("Age distance" "(absolute value)") xtitle("Expertise overlap") n(20) graphregion(color(white) margin(0)) 
graph export "${Gpath}\figures\\age_distance_overlap.pdf", replace 

* Controlling for demographic concordance index 
estimates drop _all 
local h=0 

qui reghdfe  m_30d  share_dcbo_* $XVar0 $Xvar0_phy overlap if z_demo_similarity !=., a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsing
eststo M`h++' 

qui reghdfe  m_30d  share_dcbo_* $XVar0 $Xvar0_phy z_demo_similarity  , a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsing
eststo M`h++' 

qui reghdfe  m_30d  share_dcbo_* $XVar0 $Xvar0_phy overlap z_demo_similarity  , a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsing
eststo M`h++' 

qui reghdfe  m_30d  share_dcbo_* $XVar0 $Xvar0_phy overlap z_demo_similarity c.overlap#c.z_demo_similarity  , ///
a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsing
eststo M`h++' 

esttab M* using "${Gpath}\\table\\controlling_demographic_concordance.csv",   ///
keep(overlap  z_* c.overlap#c.z_demo_similarity) replace star(* 0.1 ** 0.05 *** 0.01) se(5) b(5) staraux brackets   ///
compress   begin(";")  delimiter(";") end(";")  scalar(Mean N )


********************************************************************************
* Specialty vs Knowledge overlap 
********************************************************************************
use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1
zscore proc_overlap overlap  mean_overlap mainphy_overlap medi_overlap  overlap_cosine 

binscatter z_proc_overlap z_overlap, xtitle("Overlap in specialties" "(z-score)") ///
ytitle("Overlap in procedure proficiency" "(z-score)") graphregion(color(white) margin(0)) n(20)

graph export "${Gpath}\figures\\knowledge_specialty_overlap.pdf", replace 

*-- histogram 
reghdfe proc_overlap share_dcbo_*, a(cpfxcnesxyear) keepsingletons res(res_z)

hist res_z, frac graphregion(color(white) margin(0)) xtitle("Procedure overlap" "(residualized)")
graph export "${Gpath}\figures\\procedure_overlap_hist.pdf", replace 


*****************************************************************************************************************
* Permutation tests
*****************************************************************************************************************
use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1

sort sp_cnes procedUnic cpfxcnesxyear
* For X  
set seed 12345 
forvalues i=1/500 {
    bysort cpfxcnesxyear : gen orign = _n
    gen double r1 = runiform()
    gen double r2 = runiform()
    sort cpfxcnesxyear   r1 r2
	bysort cpfxcnesxyear  : gen p_X`i'=overlap[orign]
	cap drop r1 r2 orign
}


* For y -- reverse the seed 
set seed 54321 
forvalues i=1/500 {
    bysort cpfxcnesxyear: gen orign = _n
    gen double r1 = runiform()
    gen double r2 = runiform()
    sort cpfxcnesxyear r1 r2
	bysort cpfxcnesxyear: gen p_Y`i'=m_30d[orign]
	cap drop r1 r2 orign
}


global F0 "cpfxcnesxyear hospitalxdow hospitalxmonth"

keep p_Y* p_X* overlap m_30d   share_dcbo_* $XVar0 $Xvar0_phy ${F0} 
hdfe p_Y* p_X* overlap m_30d   share_dcbo_* $XVar0 $Xvar0_phy , a( ${F0} ) keepsingletons gen(res0_)

********************************************************************************
* Permutation X 
********************************************************************************
set matsize 5000
cap mat drop betamat
cap drop p_value 
cap drop nullbeta

forvalues i=1/500 {
dis as red "Placebo `i'"

qui reg res0_m_30d res0_share_dcbo_* $XVar0_res0 $Xvar0_phy_res0 res0_p_X`i'
	cap confirm matrix betamat
	if _rc!=0 {
	matrix betamat = _b[res0_p_X`i']
	}
	else if _rc==0 {
	matrix betamat = betamat\ _b[res0_p_X`i']
	}
}

matrix colnames betamat = nullbeta
svmat betamat, names(col)

qui reg res0_m_30d res0_share_dcbo_* $XVar0_res0 $Xvar0_phy_res0 res0_overlap
local realbeta=_b[res0_overlap]

qui: summ nullbeta
loc xmin = min(r(min), `realbeta')
loc xmax = max(r(max), `realbeta')
cap drop bin 
gen bin = round(100 * (nullbeta - r(min)) / (r(max) - r(min)))
cap drop modebin
egen modebin = mode(bin), maxmode
qui: count if bin == modebin & bin !=.
loc ymax = ceil(100 * (r(N) / 1000))

* get p value
cap drop p_value
gen p_value=abs(nullbeta) >=abs(`realbeta') 
replace p_value=. if  nullbeta==.
qui sum p_value
loc pval: di %04.3f =r(mean)
di "This is the p value: `pval'"
local pval = "`pval'"

* plot histogram
loc xlim = 1.4* max(abs(`xmin'), abs(`xmax'))
loc x_text = `realbeta'*1.25
loc y_text = `ymax' 
twoway (histogram nullbeta, fcolor(gs9) lcolor(black) percent bin(100)) , xscale(r(-`xlim' `xlim')) xlabel(#5) ///
graphregion(color(white)) xtitle(Coefficient) legend(off) text(3 `x_text' "p=`pval'", color(red) size(large)) ///
xline(`realbeta', lpattern(solid) lc(red) lw(medium)) title("Panel A. Permuting the expertise overlap", size(medium) pos(11))
graph export "${Gpath}\figures\\permutation_overlap.pdf", replace


********************************************************************************
* Permutation Y 
********************************************************************************
cap mat drop betamat
cap drop p_value 
cap drop nullbeta

forvalues i=1/500 {
dis as red "Placebo `i'"

qui reg res0_p_Y`i' res0_share_dcbo_* $XVar0_res0 $Xvar0_phy_res0 res0_overlap
	cap confirm matrix betamat
	if _rc!=0 {
	matrix betamat = _b[res0_overlap]
	}
	else if _rc==0 {
	matrix betamat = betamat\ _b[res0_overlap]
	}
}

matrix colnames betamat = nullbeta
svmat betamat, names(col)

qui reg res0_m_30d res0_share_dcbo_* $XVar0_res0 $Xvar0_phy_res0 res0_overlap
local realbeta=_b[res0_overlap]

qui: summ nullbeta
loc xmin = min(r(min), `realbeta')
loc xmax = max(r(max), `realbeta')
cap drop bin 
gen bin = round(100 * (nullbeta - r(min)) / (r(max) - r(min)))
cap drop modebin
egen modebin = mode(bin), maxmode
qui: count if bin == modebin & bin !=.
loc ymax = ceil(100 * (r(N) / 1000))

* get p value
cap drop p_value
gen p_value=abs(nullbeta) >=abs(`realbeta') 
replace p_value=. if  nullbeta==.
qui sum p_value
loc pval: di %04.3f =r(mean)
di "This is the p value: `pval'"
local pval = "`pval'"

* plot histogram
loc xlim = 1.4* max(abs(`xmin'), abs(`xmax'))
loc x_text = `realbeta'*1.25
loc y_text = `ymax' 
twoway (histogram nullbeta, fcolor(gs9) lcolor(black) percent bin(100)) , xscale(r(-`xlim' `xlim')) xlabel(#5) graphregion(color(white)) ///
xtitle(Coefficient) legend(off) text(3 `x_text' "p=`pval'", color(red) size(large)) xline(`realbeta', lpattern(solid) lc(red) ///
lw(medium)) title("Panel B. Permuting the 30-day mortality", size(medium) pos(11))
graph export "${Gpath}\figures\\permutation_mortality.pdf", replace 

*********************************************************************************************************
* Excluding large hospitals: one to many  
*********************************************************************************************************

use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1

preserve 
bysort sp_cnes: gen hospital_size=_N 
keep hospital_size sp_cnes
bysort sp_cnes: keep if _n==1
gsort -hospital_size 
gen id_hospi_sorted=_n

tempfile A 

save `A'.dta, replace 
restore 

merge m:1 sp_cnes using `A'.dta, nogen  


gen beta=. 
gen beta_min=. 
gen beta_max=. 

forvalues i=1/15 {  
qui reghdfe  m_30d  share_dcbo_* $XVar0 $Xvar0_phy overlap if id_hospi_sorted>`i' , a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsing

replace beta=_b[overlap] if id_hospi_sorted==`i'
replace beta_min=beta-_se[overlap]*1.95 if id_hospi_sorted==`i'
replace beta_max=beta+_se[overlap]*1.95 if id_hospi_sorted==`i'
}

local N=_N 
keep beta beta_min beta_max hospital_size id_hospi_sorted
collapse (mean) beta beta_min beta_max hospital_size, by(id_hospi_sorted)

qui sum beta_min
local X=r(min)
ren id_hospi_sorted id 
ren hospital_size N 
sort id 
gen X=sum(N)
gen obs=1.2*`X'
replace N=`N'-X 


keep if id<=15
twoway 	scatter beta id , msymbol(Oh) msize(vthin) mlwidth(vthin ) mcolor(black) ///
	||	scatter obs id, mlab(N) mlabsize(vsmall) mlabangle(80) mcolor(white) ///
	||	rcap beta_min beta_max id, lwidth(vthin) lcolor(black) ///
	xlabel(,valuelabel angle(45)) ///
	xlabel( 1(1)15 ) ylabel(-0.08(0.02)0) ///
	xtitle("Hospital excluded") ///	
	legend(off) bgcolor(white) graphregion(color(white))
	
graph export "${Gpath}\figures\\excl_hosp_one_to_many.pdf", replace 	


*********************************************************************************************************
* exclusion of hospitals one by one 	
*********************************************************************************************************

use  "${Gpath}data\\estimation_sample.dta", clear 
keep if main_sample==1

preserve 
bysort sp_cnes: gen hospital_size=_N 
keep hospital_size sp_cnes
bysort sp_cnes: keep if _n==1
gsort -hospital_size 
gen id_hospi_sorted=_n

tempfile A 

save `A'.dta, replace 
restore 

cap drop id_hospi_sorted hospital_size beta beta_max beta_min
merge m:1 sp_cnes using `A'.dta, nogen  

gen beta=. 
gen beta_min=. 
gen beta_max=. 

qui sum id_hospi_sorted
local max=r(max)

forvalues i=1/`max' {  
qui reghdfe  m_30d  share_dcbo_* $XVar0 $Xvar0_phy overlap if id_hospi_sorted!=`i' ///
, a(cpfxcnesxyear  hospitalxdow hospitalxmonth) cluster(sp_cnes) keepsing

replace beta=_b[overlap] if id_hospi_sorted==`i'
replace beta_min=beta-_se[overlap]*1.95 if id_hospi_sorted==`i'
replace beta_max=beta+_se[overlap]*1.95 if id_hospi_sorted==`i'

}
 
keep beta beta_min beta_max hospital_size id_hospi_sorted
collapse (mean) beta beta_min beta_max hospital_size, by(id_hospi_sorted)

qui sum id_hospi_sorted
local X=r(max)

ren id_hospi_sorted id 
sort id 

twoway 	scatter beta id , msymbol(Oh) msize(vthin) mlwidth(vthin ) mcolor(black) ///
	||	rcap beta_min beta_max id, lwidth(vthin) lcolor(black) ///
	xlabel(,valuelabel angle(45)) ///
	xlabel( 0(10)`X' ) ylabel(-0.08(0.02)0) ///
	xtitle("Hospital excluded") ///	
	legend(off) bgcolor(white) graphregion(color(white))	
graph export "${Gpath}\figures\\excl_hosp_one_by_one.pdf", replace 	

